Equilibration and thermalization of the dissipative quantum harmonic oscillator 

in a non-thermal environment 
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Institut fur Physik, Ernst-Moritz-Arndt-Universitdt, 17487 Greifswald, Germany 

We study the dissipative quantum harmonic oscillator with general non-thermal preparations of 
the harmonic oscillator bath. The focus is on equilibration of the oscillator in the long-time limit 
and the additional requirements for thermalization. Our study is based on the exact solution of 
the microscopic model obtained by means of operator equations of motion, which provides us with 
the time evolution of the central oscillator density matrix in terms of the propagating function. 
We find a hierarchy of conditions for thermalization, together with the relation of the asymptotic 
temperature to the energy distribution in the initial bath state. We discuss the presence and absence 
of equilibration for the example of an inhomogeneous chain of harmonic oscillators, and illustrate the 
general findings about thermalization for the non-thermal environment that results from a quench. 

PACS numbers: 05.30.-d 



I. INTRODUCTION 

Equilibration can be defined as the evolution of a sys- 
tem out of equilibrium towards a stationary state in 
the long-time limit. For quantum systems, the question 
arises how equilibration is possible in spite of the linear 
and unitary time evolution, how the stationary state de- 
pends on the initial conditions, and to which extent it 
can be described as a thermal state. 

General arguments relate equilibration to dephasing of 
quantum states [H-Q • Starting from the expansion of an 
initial state |^(0)) = 2 n =i V'nl 71 ) m the eigenstates \n) of 
the Hamiltonian H = £^Lj E n \n)(n\, the time evolution 
of an operator expectation value is given by 



(A(t)) = mum)) 



N 



r m ^nS E -- E ^{m\A\n) 



(1) 



m,n— 1 



In the thermodynamic limit JV->oowe can expect that 
only diagonal terms m — n survive for t — > oo, such that 
the long-time limit of the expectation value is 



lim (A(t)) 

t— ^±oo 



tr 



\ Poo A] (N ->• oo) , 



(2) 



with the density matrix = Yln=i IV'nl 2 ! 71 ) ( n \- This 
argument can be justified with the Riemann-Lebesgue 
lemma @ that states 



lim 

t— ^±oo 



f{u) e^'dw = 



(3) 



for any intcgrablc function /(w) (here: the density of 
states D(uj)). Although this argument explains the origin 
of equilibration, not much is learned about the properties 
of the stationary state p^. Especially the question of 
thermalization is left open. 
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In this paper we study equilibration and thermaliza- 
tion of dissipative quantum harmonic oscillators, using 
the standard model of a central oscillator coupled to a 
harmonic oscillator bath. For this example we can de- 
termine the stationary state poo explicitly and analyze 
its dependence on the initial conditions completely. Cru- 
cially, we allow for arbitrary non-thermal bath prepara- 
tions in our study. Thermalization is subject to addi- 
tional conditions in this more general situation, and we 
show how the temperature of the asymptotic stationary 
state is obtained from the initial energy distribution of 
the oscillator bath rather than from the initial bath tem- 
perature. We also include the case study of an interaction 
quench in an infinite harmonic chain, where undamped 
oscillations can prevent equilibration at strong damping. 

The dissipative quantum harmonic oscillator is stud- 
ied extensively in the literature I7H10| . covering such di- 
verse topics as Brownian motion [Ill4l5| . quantum fluc- 
tuations PH, driven dissipative systems (l7jh entangle- 
ment [l8[ , the existence of local temperatures [19| , or the 
second law of thermodynamics [20| . Reviews are given, 
e.g., in 21-23]. With an exact solution this model is 
also an important example for the derivation of master 
equations [24l - [27j , the discussion of fundamental statisti- 
cal relations such as fluctuation-dissipation theorems [28[ 
and their connection to detailed balance and the Kubo- 
Martin-Schwinger condition ]2£§, or for the assessment 
of numerical methods that provide a perspective for non- 
linear models [13, [H|- It appears, however, that the ques- 
tions addressed here have not been previously analyzed in 
detail, especially not for non-thermal bath preparations. 

To obtain our results we proceed as follows. After 
introduction of the model in Sec. [TTl we construct the 
exact solution for non-thermal initial states in Sec. Mil 
including the propagating function in Sec. 1111 Dl Further 
details, including the extension to driven oscillators, are 
given in App. [Aland App.|B] The central results for equi- 
libration and thermalization are formulated in Sec. IIVI 
We discuss these results for the example of an infinite 
chain of harmonic oscillators in Sec. |Vl before we con- 
clude in Sec. IVTl 
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II. THE MODEL 

The Hamiltonian for the dissipative quantum harmonic 
oscillator, 

H = H S + H B + H SB , (4) 
is the sum of the contribution of the central oscillator, 
1 



H s = 



the contribution of the harmonic oscillator bath, 



1 



H R = - 



N 



2 



\P 2 V + ulQl 



and the linear interaction term 

N 

Hsb = Q ^ XvQv 



(5) 



(6) 



(7) 



In these expressions, Q u , P v are position and momen- 
tum operators with canonical commutation relations, e.g. 
[Q/j,, P v \ = iSfj, v . Summations over Greek indices, used for 
bath oscillator operators Q v , P„, run from 1, . . . , N. We 
suppress an index for the central oscillator operators. 

The size of the coupling constants \ v is restricted by 
the positivity condition 



JY 



ft 



\ 2 



(8) 



It guarantees that the normal modes of the total 
Hamiltonian H have real frequencies, such that H is 
bounded from below [llj. A positive Hamiltonian 
can always be obtained through addition of the term 

(i/2)Ef=i( A «/^) 2( 9 2 > which leads t0 

renormalization 

of the central oscillator frequency [22(. We prefer the 
present form of the Hamiltonian since it allows for a more 
natural treatment of the harmonic chain in Sec. [V] 

Of primary interest to us is the central oscillator den- 
sity matrix 



Ps(t) =tr s [exp(-Lfft)p(0)exp(Lfft)] , 



(9) 



which is obtained from the initial state p(0) through 
propagation with the total Hamiltonian H and subse- 
quent evaluation of the partial trace tr# over the bath 
degrees of freedom. A natural choice for p(0) are factor- 
izing initial conditions 



P (0) = Ps(0)®Pb(0) , 



(10) 



which correspond to the picture that at t = the pre- 
viously isolated central oscillator is brought into contact 
with the oscillator bath. 

The restriction to factorizing initial conditions is not 
essential for the following derivations, especially not for 
the long-time limit in Sec. lIVi but it is a natural assump- 
tion that simplifies the presentation. For example, mixed 
central/bath oscillator terms drop out of the expressions 
for the central oscillator variance (see Sec. MI B|) . 



III. SOLUTION OF THE DISSIPATIVE 
QUANTUM OSCILLATOR FOR GENERAL 
INITIAL CONDITIONS 

The central oscillator density matrix ps(t) can be ob- 
tained in various ways, e.g. through transformation of H 
to normal modes @,[ll| or by using path integrals pHI^ 
based on the Feynman- Vernon influence functional for- 
malism [33-35]. The arguably simplest approach is the 
direct solution of the Heisenberg equations of motion for 
the operators Q(t), P(t), which reduces to the solution 
of a classical equation of motion. The initial conditions 
ps(0) and ps(0) enter only the evaluation of central os- 
cillator expectation values, such that we can allow for 
general initial bath states. The full solution is then given 
by the propagating function. 



A. Reduction to classical equation of motion 

As further detailed in App.fAl the central piece of infor- 
mation is the solution u(t) £ R of the classical equation 
of motion 

il{t) = -n 2 u(t)+ [ K(t-T)u(r)dr , (11) 
Jo 

which is subject to the conditions 

1. u(t) solves Eq. (|njl for t > 0, 

2. u(t) = for t < 0, 

3. the initial conditions are u(0) — 0, u(0) = 1. 
We here introduced the damping kernel 



A? 



Kit) = V — sinuU 



(12) 



The function u(t) can be calculated as the Fourier trans- 
form n 

2 f°° 

u(t) = - sinwtImF(a; + iO + )da; (13) 

7T Jo 



of the function 



F(z)= ft 2 -z 2 + £^i- 



(14) 



writing F(w + i0 + ) = lim^o^o F(oJ + irj). We note that 
the positivity condition (jHJ implies that the poles of F(z) 
occur on the real axis, such that u(t) is a quasiperiodic 
function for finite N while u(t) — > for t — > oo is possible 
in the thermodynamic limit N — > oo. An explicit exam- 
ple for the computation of u(t) is given for the harmonic 
chain in Sec. [V] (see Eq. ([79])). 
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To proceed, we introduce the partial Fourier trans- 
forms 



u(t,u) = e itul / u(t) e-' 1UT dr , (15) 
Jo 



v(t,uj)=e ltuj ii{T)e- iujT dr = u(t) + iam(t, w) , (16) 
Jo 



and define the matrices 



Tim - W 



U(t,w) 



Reu 



CO 



lmv(t,oj) 



(18) 



7 



We now obtain the central oscillator operators from the 
matrix equation 



(?§)= u(f) (*)-|:^-»(S 



Q„(0) 

(0) 



(19) 



B. Central oscillator expectation values 



Eq. (|19l) gives the operators P(i) as linear com- 

binations of the operators Q(0),P(0) and Q„(0), P„(0). 
This allows us to express central oscillator expectation 
values for t > in terms of the initial expectation values 
at t = 0. 

The linear expectation values are given by the matrix 
equation 

X(t)= (j^) = U(t)X(0)+I(t), (20) 



with the same shape as Eq. (|T9|) . In addition to the initial 
expectation values X(0) it contains the contribution 



I(t) 



I P (t) 



JV 



5^A w U(t,w w )X„, (21) 



where we mark the initial bath expectation values 

x„ : ( \i i (22) 



(Q,(o)) 
<p,(o)) 



with a breve as a notational convention. Note that if 
Xj, = 0, e.g. for a thermal bath, the 'noise term' I(i) 
vanishes. Then, position (Q(t)) and momentum (P(t)) 
of the central oscillator follow the classical equation of 
motion (|11[) . 



For the quadratic expectation values we define the vari- 
ance of operators A, B as 

Z AB = \{AB + BA) - (A) (B) , (23) 

which simplifies to Yjaa — (A 2 ) — (A) 2 for A — B, and 
write T,ab(£) — ^A(t)B(t)- We combine the central oscil- 
lator variances into the real symmetric matrix 



S(i) 



S QQ (t) s Q p(t) N 



^Sgp(i) Spp(t) y 
and denote the initial bath variances with the matrix 



(24) 



'£q„ Qm (0) E Q „p,(0) N 
Eq^(0) Ep„p m (0). 



(25) 



Note the index swap in the off-diagonal elements, and 
recall that mixed central oscillator/bath variances such 
as Eqq^ vanish for our choice (|10[) of factorizing initial 
conditions. 

We now obtain with Eq. (|19p the matrix equation 



(26) 



Similar to Eq. (|20[) . the first term results from the time 
evolution of the central oscillator according to the clas- 
sical equation of motion (jlip . and appears in the same 
form for an isolated oscillator. Only the second term 



£(t) = U(t)£(0)U T (t) + C(t) . 



C(t) 



CQp{i) Cpp(t) / 

A I ,A M U(i,w,)S^U T (t,^) (27) 



N 

E 



depends on the initial bath oscillator variances E^. 
Mixed terms in U(i), U(t, a;,/) do not appear for fac- 
torizing initial conditions. 



C. The thermodynamic limit 

Because u(t) is a quasi-periodic function for a finite 
number TV of bath oscillators, equilibration becomes pos- 
sible only in the thermodynamic limit N — > oo. We as- 
sume that for N — > oo the density of states 



1 N 



(28) 



converges to a continuous function. Note that D(cj) = 
for uj < since the bath oscillator frequencies are posi- 
tive. The coupling constants appear in the damping ker- 
nel K(t) and in Eq. (fTT|) as X 2 , and must thus scale as 
N^ 1 / 2 . We assume that 



(29) 



K = \(u v )/VN 
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with a continuous function A(w), and introduce the bath 
spectral function 



7(w) = D{uj) 



AM 2 



(30) 



with 7(0;) = for co < 0. The damping kernel is now 
given as 



(31) 



= / 7(w) sinwi dw , 
Jo 

and the positivity condition reads 



n 2 > / iM dw 



(32) 

/o w 

The function F(z) in Eq. ([T3]) for u(t) can be written 



F(z)=(tf-z>+ [~ J^foY 1 . (33) 

Under mild assumptions, the evaluation of the w-integral 
in this equation is possible by contour integration and 
results in 



F(z) = [0, — z + T(z) 



(34) 



for Imz > 0, where the complex function T(z) with 
7(0;) = =f(2/7t) Imr(iu + i0 + ) is the analytic contin- 
uation of 7(0;) into the upper half of the complex plane 
(see Sec. [V] for an example). For future use in Sec.lIVIwe 
note the relation 7(0;) \F(uj)\ 2 = (2/tt) ImF(w + iO + ) that 
follows from this representation. The analytic properties 
of F(z) determine the behavior of u(t) in the long-time 
limit, which is essential for equilibration (see condition 
(E0) in Sec. HV| : It is u(t) ->• for t ->• 00 if and only if 
F(z) has no isolated poles. 

The linear expectation values X„ enter Eq. (|2"Tj) with 
the prefactors A„ cx iV -1 / 2 . To obtain a finite result for 
the sum over N terms, also X y has to scale as iV -1 / 2 , 
which leads to the ansatz 



1 



X(w„) 



(35) 



with a continuous vector-valued function X(w). Then, 
Eq. (j2Tj) becomes 



£>(w)A(w)U(i,cj)XHdw . (36) 



1(f) 



The variances E„ M enter the sum in Eq. (|27p with the 
prefactors A^A^, cx iV^ 1 . We must now distinguish be- 
tween the ./V 2 off-diagonal terms v ^ fi, which require 
an additional 1/N prefactor for convergence, and the N 
diagonal terms v = p. Therefore, we make the ansatz 

£„ M = ^S (2) K, W)I ) + ±^(LJ u )6 Vfl (37) 



with continuous matrix- valued functions £( 2 )(cji, oj 2 ) 
and £W(w). Then, C(f) from Eq. ([27]) is the sum of 
the off-diagonal term 



C^(t)=JJ D(w 1 )£>(o; 2 )A(a; 1 )A(a; 2 ) 

x U(t, ux)±W (wi, w 2 )U T (f , W2) dwi dw 2 (38) 
and the diagonal term 



C«(f) 



W7(o;)U(t, w)S (1) (w)U T (t, w) do; . (39) 



If the initial bath state is uncorrelated, such as for a 
thermal bath or a general product state ps(0) = jOs(O)® 
' " " ® J°b(0)> the off-diagonal term C( 2 )(i) vanishes. 

When we construct the propagating function in the 
next subsection, we will conveniently assume that the 
initial bath state ps(0) is a Gaussian state. For the long- 
time limit, the situation of interest here, this assumption 
can be justified in the thermodynamic limit on general 
grounds [371 . |38| . The principal mechanism is illustrated 
with counting arguments of the following kind: Consider 
an uncorrelated bath state, where only N diagonal terms 
contribute in any sum over the bath oscillators. If we 
consider a higher order cumulant of bath operators, say 
Qs(v) = (Ql) ~ 3(Q 2 )(Q„) + 2(Q„) 3 as mentioned be- 
fore, it appears with a prefactor A 3 oc N~ 3 / 2 . There- 
fore, the total contribution of these cumulants scales as 
N x 7V~ 3 / 2 = N^ 1 / 2 and vanishes in the limit N -» 00. 
Similar counting arguments can be given for cumulants 
involving two or more bath oscillators in the presence of 
correlations. Because higher order cumulants vanish and 
only linear and quadratic bath expectation values sur- 
vive the N — > 00 and t — > 00 limit, we can treat the bath 
state as Gaussian in any calculation of the central oscil- 
lator density matrix. For the formulation and proof of a 
strict result, which is involved even under some simplify- 
ing assumptions, see (38j | . 



D. The propagating function 

Knowledge of the expectation values X(f), E(f) does 
not suffice to obtain the central oscillator density matrix 
ps(t), unless we restrict ourselves completely to Gaus- 
sian oscillator states (cf. Eq. (|43|) below). Otherwise, 
the general solution is given by the propagating function 
</(•) that, in position representation, expresses the den- 
sity matrix ps(q,q',t) = (q\ps(t)\q') for t > as 



Ps(qf,q'f,t) = // J(qf,q'f,qi,qlt)ps(q t ,q' l ,0)dq i dq' i . 

J J — 00 

(40) 

This expression must hold for all ps (0) and t > 0, and a 
fixed initial bath state ps(0). 

The propagating function can be calculated using path 
integrals and the result for a thermal bath is given, e.g., 



5 



in [21] . Within our approach it is more natural to con- 
struct the propagating function directly, using only that 
an initial Gaussian state of the joint central/bath oscil- 
lator system remains a Gaussian state during time evo- 
lution with the bilinear Hamiltonian H. With respect 
to the final remarks in Sec. IIII C[ we assume a Gaussian 
bath state ps(0). We can then consider the most general 
ansatz for J(-) that maps an initial Gaussian state ps(0) 
in Eq. (|40|) onto a Gaussian state ps{t) for t > 0, and will 
find that the parameters of this ansatz are fully specified 
through the linear maps ([20 ]) , ([25 ]) of X(i), 53 (i). The 
result is valid for arbitrary ps(0) in Eq. (|4U|) . but we do 
not need to consider non-Gaussian ps(t) explicitly. 

To translate this argument into equations we work with 
the Wigner function [39|, HI 



W(q,p,t) 



1 

2^ 



ps\q 



s s 



(41) 



defined by the relation 

W(x,t)= [ J w (x,x,t)VF(x,0)dx , (42) 
Jr 2 

where we write W(x,t) — W(q,p,t) with x = (q,p) T 
and dx = dq dp for abbreviation. Note that W(x, t) and 
Jiy(x, x, t) are real functions. 

A Gaussian state to given X(t), 53 (i) has the Wigner 
function 



W g {x,t) = 



exp[-i(x-X(t))-S- l (t)(x-X(t))] 



27r^detS(<) 

(43) 

and the most general expression for Jvk(') that respects 
this structure is an exponential function of the 14 lin- 
ear and quadratic terms in the coordinates q, p, q, p. The 
normalization L 2 W(x,t)dx = 1 of Wigner functions im- 
plies the condition 



/ Jw{x 1 x, t) dx = 1 
Jm 2 



(44) 



instead of the density matrix ps(q, q' ', t) in position repre- 
sentation (see also Refs. [H,|4l[ for a related calculation). 
The propagating function J\y (x, x, t) = Jw (q, p, q, p, t) is 



on the propagating function, which fixes the prefactors 
of the 5 terms q 2 ,p 2 1 qp,q,p in the initial coordinates. 
This leaves 9 free parameters that have to be fixed in 
accordance with the linear transformations ([20]) . (|26[) of 
expectation values. The final result is 



J w (x,x,t) = 



exp 



(x - U(t)x - I(i)) • C- X (t)(x - U(t)x - I(t)) 



27r- v /detC(t) 



(45) 



where the 4 + 3 + 2 = 9 parameters are the entries of 
the 2x2 matrix U(i) from Eq. ([T7| . the symmetric and 
positive definite 2x2 matrix C(t) from Eq. (|27l) . and the 
two-dimensional vector I(t) from Eq. ([2T|) . 

In order to check that this expression indeed repro- 
duces the transformations ([2TJ]) . ([2"6")h we can express the 
expectation values at t > in terms of those at t = 
through the evaluation of simple Gaussian integrals. To 
give an example, it is 

(Q(f)) - f qW(x,t)dx 
Jr 2 

= / gJT4/(x,x,t)T^(x,0)dxdx . (46) 

The integral of qJw{x, x, t) over x is a Gaussian integral 
with a linear term, and gives 

/ qJ w (x,x,t)dx = U QQ {t)q + U QP (t)p + I Q (t) . (47) 

JR 2 

The final integration over x in Eq. (|46p . which now in- 
volves the right hand side of (|4T)) . generates the initial 



expectation values (Q(0)), (P(0)). Therefore, we obtain 
the relation (Q(t)) = U QQ {t){Q{0)) + Dq P (t)(P(0)) + 
/ Q (t) = u(t)(Q(0)) +u(t){P(0)) + I Q (t) in accordance 
with Eq. (12^1) . Following this recipe, we find that the 
given expression (f4"5"j) for the propagating function Jiy(0 
reproduces the entire transformations (|2U|) . (|2^|) of the 
expectation values X(t), 53(t), as we required. 

If C(t) — > 0, we get a representation of the distribution 
(J(x-U(i)x-I(t)) from Eq. (gSJ). In particular for t = 0, 
where U(0) = 1, 1(0) = in addition to C(0) = 0, we 
have the correct result Jw{x 7 x, 0) = 5(x— x) in Eq. (j42|) . 

We note that the conveniently simple derivation of 
Jw{-) relies on the use of Wigner functions. Of course, 
the expressions for ps(qj, q'f, t) in position representation 
often reported in the literature can be recovered from 
Eq. ([45]) (see App.[Bj|. 



IV. EQUILIBRATION AND THERM ALIZATION 

The results from the previous section allow us to study 
the behavior of the central oscillator density matrix ps (t) 
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in the long-time limit t — ¥ oo. We can classify the be- 
havior according to the general criteria of equilibration 
and thermalization. Equilibration means convergence to 
a stationary state as expressed in the two conditions 

(El) the central oscillator density matrix ps(t) converges 
for t — > oo. 



(E2) the stationary state pi? 
dent of ps(0) ■ 



lim t _ ) . 00 ps{t) is indepen- 



Note that pi? will depend on the initial bath state 
Pb(0). Note further that the above definition of equi- 
libration does not distinguish between stationary equi- 
librium states and stationary non-equilibrium states with 
finite heat flows. The latter cannot arise for a single bath 
with continuous initial conditions as in Eq. ([37)) such that 
condition (El) is sufficient for the present study. 

Equilibration (El) implies convergence of central os- 
cillator expectation values for t — > oo. This, in turn, re- 
quires convergence of the matrix U(t) in Eqs. ([20]) . (1251) . 
Because the only stationary solution of the homogeneous 
differential Eq. (fTTj) is u(t) = 0, convergence of U(i) is 
equivalent to U(t) — ► or u(t) — > for t — > oo. There- 
fore, we assume in this section the condition 



(E0) u(t) -> for t 



oo 



as the prerequisite for equilibration (El). Under this as- 
sumption, we will be able to show convergence of expec- 
tation values and, building on this result, convergence of 
the central oscillator density matrix. 

In the weak damping limit, condition (E0) is equiva- 
lent to 7(f2) > (taking the thermodynamic limit for 
granted). This expresses the basic fact that equilibra- 
tion occurs through energy exchange with the environ- 
ment, which is not possible for an isolated oscillator with 
7(0) = 0. We note that a small value of 7(0) can re- 
sult in long transients that prevent equilibration over the 
observation time. 

Thermalization additionally requires that the station- 
ary state pi? is a thermal state, and we have the three 
increasingly stronger properties 

(Tl) the stationary state p^P is a thermal state, 

(T2) the stationary state is a thermal state pi? cx 
Q-Hs/Toc f fa e cen tral oscillator, 

(T3) the temperature of the stationary thermal state 
Pg 3 is independent of the central oscillator fre- 
quency. 

We will see that the stationary state is always Gaussian, 
which implies property (Tl). Property (T2) reduces to 
an equipartition condition on the central oscillator vari- 
ances that determine the Gaussian state, while property 
(T3) leads to a strong condition on the initial bath state. 



A. Expectation values in the long-time limit 

The assumption U(t) — > for t — > 00 implies that 
the terms U(t)X(0) in Eq. and U(t)£(0)U r (*) in 
Eq. (|26|) drop out of the expressions for X(t) and S(t) 
in the long-time limit. Only the terms I(t) and C(t), 
which depend exclusively on the initial bath preparation, 
can survive the t — > 00 limit: All information about the 
initial central oscillator state is lost. We can not imme- 
diately draw a conclusion about the long-time behavior 
because the functions u(t,u>), v(t,u>) from Eqs. (fl"5j) . (fT6|) 
do not converge for t — > 00. Instead, we note that u(t, ui) 
behaves asymptotically as 



Uas(t,u) 



u(t) e 



dr 



(t -> 00) . (48) 



Similarly, it follows v{t,uj) ~ iuju as (t, uj) for t — > 00 
from Eq. (fT6|). Consequently, the matrix U(t, u) behaves 
asymptotically as 



Reu as (t, uj) 



Imu as (t,uj)* 



U(t,u,)~ ( ^' ' uj I (t 00) , 

-uiJmu as (t,ui) Reu as (t,u>) / 

(49) 

and remains oscillating for t — > 00 even if u(t) — > 0. 

The contributions to the term I(t) in Eq. ([55)1 . say to 
(Q(t)), are of the form 

/•OO 

-Re/ D(uj)X(uj)u as {t,uj)XQ(uj)duj . (50) 
Jo 

The integrand depends on t through the factor e lujt from 
Eq. (|48|) . such that the integral is the Fourier transform of 
an integrable (by assumption even continuous) function 
of uj. If we recall the Riemann-Lebesgue lemma ((3]) we 
see that I(t) — > for t —> 00. Altogether, it follows that 
the position and momentum expectation values vanish in 
the long-time limit, i.e. X(f) — > for t —> 00. 

For the variances, a finite contribution can survive the 
t — ¥ 00 limit because the squares of the matrix elements 
of U(i, uj) occur in C(t). For example, the diagonal term 
C (1) (t) from Eq. contributes to £qq(£) the integral 



C Qo( t ) = UJ ~f(uj)cQQ(t,uj)duj 



(51) 



of the function 



c QQ {t,u) = [Reu(t,uj)] 2 t^ Q (uj) 



+ 



2\Reu(t,uj)]\hau(t,uj)] ^n) 



[lmu(t,uj)} 2 
—2 ^ppM 



Here we write, using the notation from Eq. (|2"4"]) , 
Cqq (*) Cqp(*)\ 



C«(t) 



v^qpW ^pp(*)> 



(52) 



(53) 
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for the matrix elements of C^'(t) and 



J QQ 

'<!) 
J QP 



(54) 



for the matrix elements of S^^w) from Eq. (|37|) . 
The contribution from the first term in cqq(4, a;) is 



W7(cj)[Re{i os (t, w)] 2 Sqq(w) dw 



(55) 



If we expand the square [Re u as (t, ^)] 2 according to 

[Ree ILJt z] 2 = + r ' cos 2ut - z r z t sin 2ut , (56) 

for a complex number z with z r = Re z, Zi = Im z, we see 
that in the above integral a contribution |w as (i, w)| 2 /2 
remains finite for t — > oo, while the oscillatory terms 
with cos 2u;i, sin 2uit vanish according to the Riemann- 
Lebesgue lemma Similar expressions are obtained 
for the remaining terms in C^(t). 

The off-diagonal term C^(t) from Eq. ([55)1 is given 
by a double Fourier integral and contains only oscillatory 
terms in the two frequencies u>2- Therefore, C' 2 ' (t) — > 
for t — > oo. 

We can now collect the finite contributions from the 
different terms in C^(t), to find that the central os- 
cillator variances converge to stationary values S°° = 
lim t _ ! . 00 £(i) in the long-time limit. They are given by 



y-*oo 



n PP — 



7(w) 



e ir "u(T) dr 



e 1TU1 u(T) dr 



(57) 



uiE(u) dtu , (58) 



where 



E{u) = \^t% Q {Lo) + t%{u) 



(59) 



(60) 



Comparison with Eqs. (|13p . gives the alternative ex- 
pressions 

S QQ^/°° Im ^ + i0+ )^ dw ' ( 61 ) 



Zjpp — 



ImF(w + iO + )w£(w)dw. (62) 



Recall that F(u) + i0 + ) is a continuous function according 
to our assumption u(t) — > 0. 



As noted before, the values S°° are independent of 
the initial central oscillator state. Furthermore, the ini- 
tial bath state pb(0) occurs only through the frequency- 
resolved energy distribution £{oS). In particular, the 
known equations for thermal baths Q are recovered 
whenever £{uj) = E(T,lu), where 



E(T, ft) = - coth — 
y 1 ' 2 2T 



(63) 



is the energy of a thermal oscillator at temperature T. 
Because there are no separate conditions on the two func- 
tions Sqq(oj), Epp(cj), thermalization is possible also in 
non-thermal environments (see below). 

Eqs. (|57)) - ([59|) follow directly if we assume a thermal 
bath from the outset, with initial conditions cj 2 £qq(cj) = 

tp%(u)) = E(T,u) and £% P (u>) = 0. Equipartition of 
energy allows us to combine the terms in Eq. (|52"|) to 
CQQ(t,u>) = \u(t,u>)\ 2 E(T,uj)/uj 2 , which depends only on 
the modulus of u(t,cu). We can then drop the exponen- 
tial factor e lut from Eq. (|4"51) , and convergence of £(i) 
is evident. This short cut is not available in the general 
case. 



B. Equilibration of the central oscillator 

If the initial bath state /?b(0) and the central oscil- 
lator state ps(0) are both Gaussian states, the central 
oscillator density matrix ps(t) is Gaussian for all t > 0. 
Then, ps(t) is completely determined by the values of 
X(t), E(i), and their convergence suffices to establish 
equilibration (El), and also (E2), in this case. 

Otherwise, for non-Gaussian initial states ps(0), we 
can use the propagating function Jyy(x,x, t) from 
Eq. (|45)) to find ps(t) for t — > oo. Recall that according 
to Sec. IIII O we can assume that the initial bath state is 
Gaussian in the thermodynamic limit, which allows for 
the construction given in Sec. IIII Dl 

Equilibration follows now from the observation that 
Jv/(x, x, t) converges for t — > oo whenever X(i), £(£) 
converge. The long-time limit 



lim Jv^(x, x, t) 



exp 



-a- (s c 

2 V 



2 7 rv / dct S c 



(64) 
I(t) = 
. Because 



is obtained through substitution of limt 
and limt^oc, C(t) = £°° from Eqs. ([57])— ([5" 
U(i) — > 0, the result does not depend on x. 

The long-time limit of the Wigner function Wg°(x) = 
lim^oo Ws(x, t) follows immediately with Eq. pSl) : The 
integration over x in the resulting expression 



Wf (x) = / (x) W s (x, 0) dx = (x) 

JR 2 



(65) 



evaluates to one because VFs(x, 0) is normalized, such 
that Wg°(x) is equal to J^(x). In other words, the sta- 
tionary state pi? is a Gaussian state (|4"3")l with parameters 
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X = 0, S = E°°. These parameters depend on the ini- 
tial bath state according to Eqs. ([57)1 , (|55|) . but they are 
independent from the initial central oscillator state. This 
proves equilibration (El) and (E2) for general initial cen- 
tral oscillator states. In particular, the stationary state 
is Gaussian also for non-Gaussian initial states. 

We note that the propagating function in position rep- 
resentation does not converge in the long-time limit (cf. 
App.[B]), which prevents an equally simple argument. 



C. Thermalization of the central oscillator 

Because the stationary state p?? in the long-time limit 
is a Gaussian state for which only £qq, Ep°p are non- 
zero, it can always be interpreted as the thermal equilib- 
rium state of some harmonic oscillator. This establishes 
the weakest thermalization property (Tl). 

The effective oscillator frequency £1^ and temperature 
Too associated with pf are 



the condition 



fit = 



J pp 



T — 

-L DO 



-arcoth 



n QQ n PP 



(66) 

Generally, Ooo is not equal to the central oscillator fre- 
quency such that the stronger property (T2) is not 
fulfilled. By Eq. (|66|) . the condition £1^ = fl is equiv- 
alent to cquipartition of kinetic and potential energy 
(P 2 ) = Epp = n 2 Z^ Q = fl 2 (Q 2 ). The violation of 
this condition arises from the integrations over w in 
Eqs. ([57]), ([58]) or ([61]), ([62]), which cover a finite energy 
range and include values lo ^ Q. Quantum corrections of 
this type are characteristic for strong damping Q. 

Equipartition of energy is achieved in the limit of weak 
damping (7(0) — > 0, when according to Eq. (IM)) the 
function (2/tt) lmF(u + i0+) in Eqs. (HI]), converges 
to 2(5(w 2 - n 2 ) = (6(uj + 8(u - H))/fi. Therefore, 
the values fi 2 £gg = £p P = £{Vi) are obtained. This 
establishes the stronger thermalization property (T2) in 
the weak damping limit. 

For these values of EgU, E^p it is (cf. Eq. 



t (wd )(0) 



-arcoth 1 X , (67 

2 n K ' 



such that the stationary state p'c? is a thermal equilibrium 
state of the central oscillator. The temperature T oc (ri) is 
determined by the energy £(f2) of the bath oscillator at 
frequency f2 in the initial state. Note that the assumption 
(E0) implies 7(f2) ^ and D(lu) ^ 0, such that the value 
of £(f2) is defined. In particular, £ (Q) > f2/2 and the 
argument of arcoth(-) is equal to or greater than one. 



Still, the asymptotic temperature T a 



-,(WD) 



(Q) 

from Eq. (|67[) is a function of ft. The functional depen- 
dence is determined by the choice of £{oS). If we demand, 
for the strongest thermalization property (T3), that 
is independent of fl we have to solve Eq. (|r?T|) to obtain 



f» = |coth^- 



(68) 



Note that this is a condition on the particular combi- 



nation £{uj) of initial bath variances Sgg(w), Epp(w), 
and not on the individual functions. Therefore, any ini- 
tial bath preparation with £{uj) = E(Tq,ui) results in the 
same stationary states as the thermal bath at tempera- 
ture To. One example for this additional freedom is the 
choice 



^(1) 



_ coth(^/2r ) - 1/2 
^qqM 



t<p P (u) = w/2 , 
(69) 

and arbitrary sL^w). It can be realized, e.g., by su- 
perposition of coherent oscillator states at different po- 
sitions. This initial bath state is not a thermal state 
for To > 0, in particular it violates equipartition of en- 
ergy l^E^qH = Ep P (o;). But since £{uj) = E{T ,uj) 
we find that the stationary central oscillator state p"c? is 
identical to that obtained with a thermal bath at temper- 
ature To: Thermalization is well possible in non-thermal 
environments, even those far from thermal equilibrium. 



D. Summary 

In summary, we have a hierarchy of conditions for equi- 
libration and thermalization: 

(El), (E2) the central oscillator equilibrates whenever 
u(t) — ¥ for t — > 00, 

(Tl) the stationary state is always a Gaussian and 
thermal state, 

(T2) equipartition of kinetic and potential energy 
occurs precisely at weak damping, 

(T3) the asymptotic temperature is indepen- 
dent of the central oscillator frequency under 
the additional condition (|68l) on £{oj). 

It is a special feature of linear systems such as the 
one studied here that equilibration depends only on the 
asymptotic behavior of the solution u(t) of a classical 
equation of motion (TTTI) . Another feature is that the sta- 
tionary state always is Gaussian such that equilibration 
implies thermalization, albeit only in the weak sense of 
property (Tl). We noted earlier that in the situation 
studied here, with coupling to a single bath, a station- 
ary state does not admit finite heat flows as would be- 
come possible for several baths with different prepara- 
tions £(ui). Therefore, conditions (El), (E2) capture the 
standard notion of thermodynamic equilibrium. 

We note that a consistent definition of thermalization 
requires the strong property (T3). Suppose we deal with 
two central oscillators with frequencies Q\ 7^ f^- In the 
weak damping limit, the stationary state is the product 




• • WW • WW • WA • WW • - 

ii = -2 -1 1 2 

FIG. 1. (Color online) Sketch of the infinite harmonic chain 
as defined in Eqs. (J70]), ffTl . 



state of two independent thermal states with respective 
temperatures T OQ (f2i) and T oo (0,2). Such a state is only 
a thermal state of the combined system comprising the 
two oscillators if T"oo(^i) = ^00(^2)- Therefore, thermal- 
ization of multiple oscillators, already in the weak sense 
(Tl), requires the strong property (T3) and thus condi- 
tion (1551) (but recall that this condition can be fulfilled 
also for non-thermal environments as in Eq. (I69[l ). 



V. THE INFINITE HARMONIC CHAIN 

As an example for equilibration in a non-thermal en- 
vironment we consider an infinite chain of harmonic os- 
cillators (see Fig.[T|). Oscillators in the right (n > 1) and 
left (n < —1) half of the chain, with frequency Of, , are 
coupled to their neighbors (nil) with spring constant 
fcb. They form the harmonic oscillator bath for the cen- 
tral oscillator at n = 0, with oscillator frequency VL and 
coupling k to the oscillators at n = ±1. For ft — fij, and 
k = kf, we obtain a homogeneous, translational invariant 
chain. 

Related examples have been studied in numerous pub- 
lications, see e.g. [l3l |42T - I48j . The behavior for thermal 
initial conditions, e.g. in a homogeneous chain [42| or a 
chain with a single heavy mass [13| . is well understood. 
Equilibration in a harmonic chain with non-thermal ini- 
tial conditions as discussed in Refs. j4j| HI] can be ex- 
pressed in terms of our conditions from Sec. IIVD1 Gen- 
eral arguments for the appearance of Gaussian states in 
the long-time limit are given in [38|,|46j]. Still, a satisfac- 
tory and explicit analysis of equilibration and thermal- 
ization of the simple chain in non-thermal environments 
is missing. Some studies assume too quickly that equili- 
bration implies thermalization, in the sense of our con- 
dition (Tl), failing to note, e.g., that the appearance of 
Gaussian states is the general behavior of linear systems 
and unrelated to thermalization as expressed by condi- 
tion (T3). According definitions of 'temperature' have to 
be taken with care. In addition we must carefully analyze 
the role of undamped oscillatory behavior that prevents 
equilibration and, therefore, thermalization. 



Hr 



- OO 



n— 1 ri— 1 



h ^ qnq 



n+l 



(70) 



for the harmonic oscillator bath (with operators q n , p n 
for the oscillator at site n ^ 0) to normal modes. The 
same transformation has to be applied to the operator 
kqi in the coupling term 



H. 



SB 



~kQ(qi 



q-i 



(71) 



between the central oscillator and the chain oscillators at 
n = ±1. It suffices to treat one of the two half- infinite 
chains explicitly, say the right chain n > 1 as in Eq. (|70p , 
and include a factor of two in 7(0;) to account for the 
left chain n < — 1. Note that in doing so we implicitly 
assume identical initial conditions for both sides of the 
chain and thus exclude the possibility of stationary non- 
equilibrium states with finite heat flow between the right 
and left half-infinite chain. 

The normal modes of Hb are the standing wave solu- 
tions f v (n) oc sin {j^j) (for a finite chain of length N), 
and after a few lines of calculation we obtain the spectral 
function 



7M = §Sv /4fc 6-(^-^) ; 

7T /W 



for \n 2 b 



< 2k b 



(72) 

00. It is j{uj) — for 
> 2kb, and we impose the positivity condition 



in the thermodynamic limit 

in? 



> 2fch > to exclude negative frequencies of the bath. 
To proceed it is convenient to introduce the dimension- 
less model parameters 



Kb 



2h 
"fa 



2k 
"fa 



n 



and to use the normalized quantities 



u> = — — , t = tflb 
lib 



Note that < K b < 1. 



u(t) — flbu(i) 



(73) 



(74) 



B. Conditions for equilibration in the harmonic 
chain 

As discussed in Sec. IIV1 equilibration depends entirely 
on the decay of the function u(t) for t — > 00, and thus 
on the absence of poles in F(z) from Eq. (f3"3")) . To obtain 
F(uj), we use the representation (jM)) with the complex 
function 



A. Mapping onto the central oscillator model 



T(z) 



Ak 2 



(75) 



To address the harmonic chain within the formalism 
from Sees. Hi]— llVl we must transform the Hamilton oper- 



where the branch cut of the root must be chosen along 
the positive real axis, and the minus (plus) sign applies 
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t 1 r 



Rer' 




i i i 

1 

FIG. 2. (Color online) Real (dashed curve) and imaginary 
(solid curve) part of F(uj + i0 + ) for Kb = 1/2 and lo > 0. For 
uj 2 = {1 — Kb, 1,1 + Kb} the function value is { — 1, — i, 1} x 
(ftfift) 2 /Kb, respectively. 

for Rez > (Rez < 0). Note that the positivity condi- 
tion ©, which can now be rewritten as ft 2 +T(iO + ) > 0, 
requires that 

^>^(l-0^). (76) 

Before we can determine the function u(t) with 
Eq. (fT5|) we must consider the possibility of isolated poles 
of F(z). According to Eq. (|34l) we have to compare the 
functions u 2 — £1 2 and Rer(u; + i0 + ) in regions where 
ImT(w + i0 + ) = 0. From the qualitative behavior of 
r(cj + i0 + ), shown in Fig. [2 we deduce that isolated 
poles of F{z) do not exist if and only if the inequalities 

1 _ *tu?L <n 2 ,<i + li^jf. (77) 

Kb K b 

are fulfilled. The first inequality excludes poles in the 
interval uj 2 < 1 — Kb, the second inequality in the interval 
uj 2 > 1 + Kfc. Another more fundamental restriction is the 
positivity condition (|76j) , which is however less restrictive 
than the present condition. 

The admissible parameter combinations for equilibra- 
tion of the harmonic chain that follow from condition (|77|) 
are depicted in Fig. [3] We note the basic restrictions 

k < K b and \l-tt 2 \<K h . (78) 

The second inequality guarantees that the central os- 
cillator frequency fl r lies within the interval Q € 
[VI — Kb, VI + Kb\ where 7(a)) > 0. If this is fulfilled, 
equilibration is always possible for sufficiently small k. 
Since k& < 1, it restricts the admissible parameters to 
the rectangle (K b ,n 2 ) G [0, 1] x [0, 2]. 

Condition ((77)) is always fulfilled for the homogeneous 
chain (and we note that k — Kb requires Q r = 1). The 
chain studied by Ullersma corresponds to parameters 
Kb = 1 and Q r — k (O 2 equals the mass ratio [i in [Uj]). 
Condition (|77|) is fulfilled if fi r < 1, i.e. only for a heavy 
mass. Both examples lie on the boundary of the admis- 
sible parameter space, with one or two of the inequalities 
in (|77[) becoming equalities. 




FIG. 3. (Color online) Diagram of the admissible parameter 
space for equilibration according to condition (|77[) . The white 
triangular region above the solid black lines is the maximal 
set of allowed parameter combinations. Outside of this region 
an isolated pole exists even in the weak damping limit k — > 0. 
For k > 0, the region of admissible parameters shrinks as 
depicted by the dashed black curves. The parameter combi- 
nations of homogeneous chains (Q r = 1) corresponds to the 
cusps Kb = k of the curves, marked with red dots. The pa- 
rameter combinations of chains with a single heavy mass [13J 
correspond to the intersections of the curves with the Kb = 1 
line at Q r = k, marked with green squares. At these points, 
condition ()77|) coincides with the positivity condition (|76[l. 



C. Dynamical evolution of the harmonic chain 

Depending on parameters, the harmonic chain features 
rich dynamical behavior. For parameter combinations 
that fulfill condition (|77|) the explicit result for u(t) from 
Eq. (Ull) reads 

u(t) = — — / sinwi \ k 2 — (1 — u) 2 ) 2 



K 2 (cj 2 -n 2 ) 2 -2k 2 (q 2 -i)(lo 2 -n 2 ) + K 4 ' 

(79) 

For parameter combinations violating condition (|77|) iso- 
lated poles of F{z) occur and additional (undamped) sine 
functions £i sin Clrf must be added to this expression. Ac- 
cording to Eq. the poles of F(z) are the solutions of 
ft 2 — il 2 + T(Qi) — 0, which gives a quadratic equation for 
the harmonic chain such that zero, one, or two (positive) 
poles are possible. 

As an example let us consider the case k — 1/2. The 
restrictions on the parameters arising from the positivity 
condition (fTrl)) and the stronger condition (|77p are sum- 
marized in Fig. |4j We now follow the dashed path in 
this figure and plot the position ^ 2 / 2 °f isolated poles 
and their total weight £ = £1 + £2 in Fig. [SJ Only for 
parameter combinations in the white unshaded area in 
Fig.|H which corresponds to the part between the dashed 
vertical lines in Fig.[5l condition (|77|) is fulfilled. Accord- 
ingly, only panel (c) in Fig. [5] (the parameter combination 
"c" in Figs. [H [5J shows a situation where u(t) — > for 
t — > 00. Otherwise, one (parameter combination "d") 
or two ( "a" , "b" ) isolated poles exist if one or both in- 
equalities from (T77|) are violated. Then, the amplitude of 
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FIG. 4. (Color online) The (fi 2 ,K(,) parameter space of the 
infinite harmonic chain for k = 1/2. The solid/dotted black 
curves give the boundary of the two regions defined by each 
of the inequalities from condition (|77[) . The central oscillator 
equilibrates for parameters in the unshaded region above the 
solid curve, where the condition is fulfilled (note that Kb < 1). 
For parameters lying between the solid and dotted curves one 
of the two inequalities is violated and a single isolated pole 
of F(z) exists. Below the dotted curve F(z) has two isolated 
poles. The dashed region to the left indicates where the posi- 
tivity condition \7Q\ is violated, but such parameters already 
violate the first inequality in (|77|) . Both conditions coincide at 
Kb = 1, flr = k. The dashed red/green lines indicate the path 
followed in the next Fig. [5] the crosses marked a — d indicate 
the parameters used in Fig. [B] 



^0.5 




FIG. 5. (Color online) Position Sli and total weight £ of iso- 
lated poles of F(z). We set k — 1/2 and change Sl r , Kb along 
the dashed path from Fig. [4] i.e. from Q r = 0, Kb — 0.2 to 
Sl r — 2, Kb = 0.8. The position of the poles is compared to the 
continuum of bath modes in the interval up" G [1 — Kb, 1 + Kb], 
filling the grey area around fif = 1 in the plot. Between 
the two vertical dashed lines at il r = 1, Kb = 1/2 (left) and 
Kb = 0.8, Sl r = 1.4875 (right) no poles exist in agreement 
with condition ([77]). For SI 2 < 0.12628, in the dashed region 
to the left, the positivity condition (|76p is violated and one 
Si 2 becomes negative. 



oscillations in u(t) remains finite in the long-time limit. 

For strong damping situations (k ~ 1) shown in Fig. [6] 
the function u(t) deviates significantly from an exponen- 
tially decaying function, even in the absence of poles 
(panel (c)). Exponential decay occurs only for weak 
damping k « k^. For — 1| <C Kf, we have 



u(t) 



sin(f2 r t ) 



exp i 



(«<!), (80) 




FIG. 6. Function u(t) for the harmonic chain with k = 1/2. 
The parameters from panels (a)-(d) correspond to the crosses 
in Figs. H [S] They are: (a) K b = 0.2, SI, 2 = 0.4 (two poles 
Hi = 0.48, Sl 2 = 1.10, £i = 0.82, £2 = 0.10), (b) K b = 0.4, 
SI 2 = 1 (two poles fii = 0.76, fi 2 = 1.20, £1 = 0.26, £2 = 
0.26), (c) K b = 0.6, SI 2 = 1 (no pole), (d) K b = 0.8, SI 2 = 1.6 
(one pole fii = 1.35, £1 = 0.50). 
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FIG. 7. (Color online) Function u(i) for the inhomogeneous 
(left panel, with Sl r — 1, Kb = 0.5, k — 0.1) and homoge- 
neous (right panel, with Sl r = 1, Kb = 0.1, k = 0.1) harmonic 
chain at weak damping. The dashed red curves indicate the 
exponential decay from Eq. (|80[) (left panel) and the asymp- 
totic decay oc of the Bessel function from Eq. lf8T|l (right 
panel) . 



as plotted in Fig. [7] (left panel). Note that the case 
7(0) = 0, with an undamped sine function in the weak 
damping limit, is excluded by the second inequality in 
Eq. 478]). 

For the homogeneous chain, with Sl r = 1 and Kb = k, 
the weak damping limit gives a different result. Since 
Kb = k, the width of the continuum of bath states shrinks 
to zero for k — > such that we do not obtain exponential 
decay of u(t). Instead, it is 



s'mt (k <§; 1, horn, chain) (81) 



with the Bessel function J (x) (cf. Refs. |42| - |44j ). Ac- 
cording to condition ([77]) isolated poles of F(z) cannot 
occur in this situation. From the asymptotic behavior 
of the Bessel function we find that here u(t) decays only 
as 2(nKt)~ 1 / 2 for t ^> 1, as shown in the right panel of 
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Fig. [JJ Exponential decay in the weak damping limit is 
only achieved if the coupling k of the central oscillator to 
the chain becomes small also in comparison to the width 
(~ Kb) of the continuum of bath states. 



D. Thermalization after a quench 

According to the previous discussion, the central os- 
cillator in the harmonic chain equilibrates precisely for 
parameter combinations that fulfill condition (j77[) . We 
now study, under these conditions, thermalization after 
a quench that generates a non-thermal environment for 
the central oscillator (cf. Eq. below). 



1. Initial conditions generated by the quench 

We imagine that for t < all oscillators are decoupled 
(k = Kb = 0) and in thermal equilibrium at temperature 
Tq. Every oscillator has the same variance 



il 2 b t qq (n) 



Spp(n) = E(T ,Q b ) , 



(82) 



and we do not need to specify further initial expectation 
values if we are only interested in the stationary state in 
the long-time limit. 

At t = we quench the system by cranking up the 
coupling to finite values k, Kb > 0. Since T, gq (n), E pp (n) 
do not depend on n, transformation to the normal modes 
of the bath results in constant functions 



(i) 
pp 



(u) = E(T ,n b ) 



(83) 



for the initial bath variances at t = 0. The initial bath 
state is uncorrelated with £( 2 )(u;i, W2) = 0. 

According to Sec. HVl the stationary state in the long- 
time limit depends only on the frequency-resolved energy 
£ (u>) of the initial bath state, which for the present ex- 
ample is given by the function 



£(w) 



1 



-E(T ,n b 



(84) 



This function acquires a dependence on uj through the 
dispersion of the bath modes after the quench, but it does 
not fulfill Eq. (|B"5|) . We thus see that the thermal equi- 
librium state of uncoupled oscillators before the quench 
corresponds to a non-thermal state of the coupled chain 
of oscillators after the quench. According to condition 
(T3) from Sec. lIVDl we expect that the temperature Too 
of the stationary state depends on the central oscillator 
frequency VL r even at weak coupling. 

From Eqs. (|57 [> -(|5g j) or Eas. ([B"T ]) . {52} , the variances 
in the long-time limit are obtained as 



Soo 
nn — 



1 

2fif 



/ 



E(T ,n b ) , 

(85) 



Vj?P = ~(l + n£)E(T ,Slb). (86) 
We will give further results using normalized quantities 
floe — floo/flb , T^ = Too/rifc , To = To/Qb ■ (87) 
choosing fib as the unit of energy. 

2. Thermalization (T2) 

We recall that according to property (Tl) the station- 
ary state is always a thermal state of some harmonic 
oscillator Hamiltonian, such that we should check the 
stronger property (T2). From Eq. (|TJ7J|) . the effective fre- 
quency associated with the stationary state is 



VL 2 



<> 2 



K 2 

TFk^ 
ll r K b 



(i->/r^?) 



We observe that equipartition of energy, i.e. fioo = £l r , 
can be achieved only in the weak damping limit k — > 0. 
For k > 0, it is always f2oo < £l r - This confirms the 
conditions given for property (T2) in Sec. IIVDI 



3. Thermalization (T3) 
For weak damping, Eqs. ([H5]l. (1561 simplify to 



fl 2 X% Q = ££p = - (1 + n 2 r ) E(T ,n b ) (for k y 0) . 

Equipartition of energy in the stationary state is evident, 
and the thermalization (T2) property fulfilled. To check 
property (T3), we calculate the temperature 



— — = arcoth 



~(fi, 



-M coth ( i 
fl r J V2T 



(90) 



of the stationary state with Eq. (|66p or the weak damping 
result (jBT)) . We see that T 00 (f2 r ) depends explicitly on the 
central oscillator frequency f2 r , as depicted in Fig. [8] It 
is Too = T only for fl r = 1. As discussed before, this 
results from the fact that £(uj) after the quench violates 
condition (f6"5|). 

We note that Kb does not appear in Eq. (|90|) . In the 
present example the value of Kb only determines the ad- 
missible values of Q r that lead to equilibration, as given 
by the second inequality in Eq. (|78p . Once equilibra- 
tion has been observed, the temperature of the station- 
ary state at weak damping depends only on the value of 
£ (f2) not on the functional dependence of the spectral 
function 7(0;). 
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FIG. 8. (Color online) Temperature Too of the station- 
ary thermal state at weak damping as given in Eq. (|90l) . 
It is shown as a function of £l r for different temperatures 
To — 0.2, . . . , 0.8 of the initial state, as indicated. Note that 
Too does not depend on Kb, but the admissible values of Q r for 
which equilibration occurs are restricted by the second con- 
dition in Eq. (|78|) (see also Fig. [3]). In particular, it must be 
< Sll < 2. 




FIG. 9. (Color online) Frequency floo (dashed curve) and tem- 
perature Too (solid curves) for the homogeneous chain, from 
Eqs. (|93[) . (|94[) and shown as a function of k. The temperature 
curves are plotted for To = 0, 0.2, 0.5, 1 as indicated. 



4- The homogeneous chain 

For the homogeneous chain with Kb = n, Q r = 1 
Eqs. (JHSJl, (ESI) simplify to 



v^OO 

2-Tin — 



1 



1 



Ep P = E(T , fib) 



E(T ,fl b ), (91) 



(92) 



Equipartition of energy is violated for any k > 0, such 
that the effective frequency 



fit, = 



l + ii-K 2 )- 1 / 2 



(93) 



associated with the stationary state deviates from the 
central oscillator frequency (it is always fioo < 1). The 



2T a 

flr^r 



= arcoth 



coth 



1 



1 + (1-K 2 ) 



-1/2 



(94) 

1/2 for k — ¥ 1 



It is Too > T for k > 0, for example T a 
and T -> (see Fig. [5). 

The situation simplifies again in the weak damping 
limit k — y 0, where we recover from Eqs. (|9T|). (|92|) the 
equilibration/thermalization result for the homogeneous 
chain formulated in Refs. [43|, [44|: At weak damping the 
central oscillator evolves into a stationary thermal state, 
with equipartition of energy f2 2 EgQ = Epp = E(T ,flb)- 
Because of translational invariance this statement applies 
to every chain oscillator. 

We note, however, that thermalization of the homoge- 
neous chain is not perfect. As discussed in Sec. lIVDl ob- 
servation of a single oscillator in the homogeneous chain 
is not sufficient to establish thermalization of the entire 
chain. Thermalization fails for a finite chain segment 
consisting of two or more oscillators, because property 
(T3) is not fulfilled as seen in Eq. ([SO]). Note that there 
is no possibility to check property (T3) directly for the 
homogeneous chain (fi r — 1 is fixed here), such that re- 
sults restricted to this situation have to be interpreted 
carefully 0,13. 



VI. CONCLUSION 

Our study of the dissipative quantum harmonic os- 
cillator addresses equilibration and thermalization in 
non-thermal environments. Equilibration is the generic 
behavior, which is prevented only in situations where 
the classical oscillator equation of motion possesses un- 
damped oscillatory solutions. The infinite harmonic 
chain is an example for this behavior. 

Thermalization of the central oscillator depends on ad- 
ditional conditions. Just as for thermal environments, 
equipartition of energy requires the weak damping limit 
but is independent of the precise initial conditions. The 
asymptotic temperature Too is obtained from the energy 
distribution £ (u>) in the initial bath state, and generally 
depends on the central oscillator frequency fi. If we de- 
mand that Too is independent of fi, another condition on 
£ (ui) follows. This condition is essential for simultane- 
ous thermalization of several oscillators, when a thermal 
state of the combined system is obtained only if the same 
asymptotic temperature is assumed by each oscillator. 

Part of the behavior discussed here generalizes to sys- 
tems with non-linear interactions. First, we note that 
equilibration is possible although the linear system is in- 
tegrable. Equilibration occurs because, in a rough sense, 
the reduced density matrix of the central oscillator in- 
volves an average over conserved quantities of the joint 
oscillator-bath system. In other words, equilibration of 
small systems embedded in a large environment does 
not require ergodicity. Second, because of the linearity 
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and unitarity of quantum mechanical time evolution the 
stationary state depends explicitly on the initial (bath) 
state. But already for the linear system some properties, 
such as equipartition of energy, are independent of the 
initial conditions. Furthermore, the stationary state de- 
pends only on the energy distribution £{uj) in the initial 
bath state. Effectively, information is lost in the long- 
time limit and thermalization is possible for a large class 
of (non-thermal) initial states. 

We did neither discuss the generalization of the 
fluctuation-dissipation relation to the present non- 
thermal setting, nor the role of stationary non- 
equilibrium states with finite heat flow that would re- 
quire coupling to at least two baths with different prepa- 
rations. Multi-time correlations functions can be com- 
puted within the present formalism, which will allow for 
the analysis of both issues in the future. 
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Inserting this result into Eq. (IA1[) gives the equation 
of motion 



Q{t) = -tt 2 {t)Q(t) + / K(t-r)Q(T)dT-N(t) (A4) 



for the central oscillator operator Q(t), with the damping 
kernel K(t) from Eq. (TT2"j) and the noise term 



N in t 

N(t) = y2x v (coau v tQ v (0) + ^^-P v (0)) . (A5) 



Eq. (|A4|) is an inhomogeneous linear integro-differen- 
tial equation, which can be solved through solution of the 
classical equation of motion 



d tt u(t,t') = -n 2 (t)u(t,t') + J K(t-T)u(T,t')dT . 

(A6) 

We need the two solutions u\(t, t'), U2(t, t') to initial con- 
ditions ui(t,t) = 1, dtu\(t,t')\t=t' = 0, and U2(t,t) = 0, 
dtU2(t, t')\t=t> = 1. The solution of the operator equation 
of motion for Q(t) is then given by 



Appendix A: Operator equations of motion and 
their solution 



The solution of the dissipative quantum harmonic os- 
cillator model through operator equations of motion in- 
stead of transformation to normal modes of H allows for 
a simple treatment of general initial conditions and time- 
dependent coefficients. We here list the relevant steps of 
the derivation omitted in Sec. Mil and allow for a time- 
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dependent central oscillator frequency Q(t) (cf. Ref. 
for a path integral calculation). 

The Heisenberg equations of motion A(t) — i[H,A(i) 



Q(t) = ui(t, 0)Q(0) + u 2 (t, 0)P(0) - / u 2 {t, t)N(t) dr , 

Jo 

(A7) 

and it is P(t) = Q(t). 

With the partial Fourier transforms 



u{t,uj) = / u 2 (t,T)e iulT dr 
Jo 



v(t,u)= / d t u 2 (t,T)e l " T dT , 
lo 



(A8) 



(A9) 



N 



Q(t) = p(t) , P(t) = -n 2 (t)Q(t) - J2 KQu{t) 

(Al) 

for the position and momentum operator of the central 
oscillator, and 

Quit) = P v (t) , Pu{t) = -co 2 Mt) - A„Q(t) (A2) 

for the bath oscillators. 

We can read Eq. (|A2[) as an inhomogeneous linear 
equation for Q u {t). Using the Green function for the 
homogeneous equation Q v (t) = —uj 2 Q v (t\ we find 

Quit) = cos u) v t Qi/(0) + — smuj v t P„(Q) 

-A„ / — sinw y (<-T)Q(r)dr. (A3) 
Jo 



and the definition of matrices 

'ui(t,0) u 2 (t,0) 



U(t) = 



v atui(*,0) d t u 2 (t,0). 



(A10) 



and XJ(t,uj) as in Eq. P^)l the operators Q(t), P(t) are 
given by the matrix Eq. (1191) . 

For constant f2(£) = f2, the function u(t) used in 
Sec. |ni]is recovered as u(t) = u 2 (t,0), and it is ii(t) = 
Ui(t, 0) (while ui(t,t') ^ d t u 2 (t,t') for time-dependent 
f2(i)). Then, the partial Fourier transforms u 2 (t,uj) and 
V 2 (t, w) are related by Eq. (IT^l) and V(t) is given by the 
simpler expression (|17p . 

Eqs. d20j), (HH and Eqs. <j26j) , (J27J) for the calculation 
of expectation values and variances and Eq. (1431) for the 
construction of the propagating function remain valid for 
time-dependent Q(t). 
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Appendix B: Propagating function in position 
representation 

The propagating function in position representation is 
the Fourier transform 

J(q f ,q' f ,qi,ql,t) = i- JJ°° e -W«-*i) 

T (Qf + v'f . q l + q' t a , 

xJw[ — ^">P! — 2"^' P '*J dpdp 

(Bl) 

of Eq. (|45|) . It results in the expression 



lie I 



2tt 



exp 



ji£ 2 +32xy + j 3 y 2 



J{Y,y,X,x,t) 

+i((j4X + j 5 y)X + {j 6 x + j 7 y)Y + j s x + j 9 yj 



(B2) 

where we write Y = (qf + q'f)/2, y = ?/ — q't, X — 
{q% + q'i)/2, x = qi ~ q[ f° r abbreviation and drop the 
time argument in jk = jk{t)- The 9 real parameters 
ji , . . . , jg in this expression are related to the parameters 



of Jw(x,x,t) in Eq. (j45|) through 



Ji 



a 



QQ 



2U 2 



h - ~t; c pp - 2 ^2 P C QQ 



J4 



1 

-( 

2 

Uqq 

Uqp 
Upp 



J2 



ui 



Cqp Upp 



U, 



QP 



U 2 

U QP 



Upp 



QP U QP 

UqqUpp 



QP 



35 = Upq - 



Uqp 



36 = 



Uqp ' 



J 8 



1 T T UPP T 

1q j 39 = ip - 77 — Iq 



Uqp 



Uqp 



(B3) 



Explicit insertion of U(t) from Eq. (fT7|) or (|A10|) gives 
expressions that allow for direct comparison with the lit- 
erature. For example, the expressions given in Ref. [l7j 
are recovered for I(t) = such that the terms j%x, j$y 
vanish. 

Obviously, the position representation leads to less 
convenient expressions for the propagating function, and 
obscures the clear formal structure of Eq. (|45j) . In par- 
ticular, the expressions (IB3|) are singular for u(t) — ¥ 
0, which gives a complicated representation of the 5- 
distribution for the propagating function at t = and 
t — > oo instead of the simple limit for Jv^(x, x, t) (cf. 
Eq. (IM])). 
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